Stochastic calculus for uncoupled continuous-time random walks 
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The continuous-time random walk (CTRW) is a pure-jump stochastic process with several appli- 
cations in physics, but also in insurance, finance and economics. A definition is given for a class of 
stochastic integrals driven by a CTRW, that includes the Ito and Stratonovich cases. An uncoupled 
CTRW with zero-mean jumps is a martingale. It is proved that, as a consequence of the martingale 
transform theorem, if the CTRW is a martingale, the Ito integral is a martingale too. It is shown 
how the definition of the stochastic integrals can be used to easily compute them by Monte Carlo 
simulation. The relations between a CTRW, its quadratic variation, its Stratonovich integral and its 
Ito integral are highlighted by numerical calculations when the jumps in space of the CTRW have 
a symmetric Levy Q-stable distribution and its waiting times have a one-parameter Mittag-Leffier 
distribution. Remarkably these distributions have fat tails and an unbounded quadratic variation. 
In the diffusive limit of vanishing scale parameters, the probability density of this kind of CTRW 
satisfies the space-time fractional diffusion equation (FDE) or more in general the fractional Fokker- 
Planck equation, that generalize the standard diffusion equation solved by the probability density 
of the Wiener process, and thus provides a phenomenologic model of anomalous diffusion. We also 
provide an analytic expression for the quadratic variation of the stochastic process described by the 
FDE, and check it by Monte Carlo. 

PACS numbers: 02.50.Ey, 05.40.Jc, 



I. INTRODUCTION 

A. The continuous-time random walk 

The continuous-time random walk (CTRW) is a pure- 
jump stochastic process used as a model for standard and 
anomalous diffusion when the sojourn time at a site is 
much greater than the time needed to jump to a new po- 
sition, i.e. when jumps can be considered instantaneous 
events. The CTRW has been introduced in physics by 
MontroU and Weiss [l|; other seminal papers on its ap- 
plication to standard and anomalous transport phenom- 
ena are due to Scher and Lax 0, Q and to MontroU and 
Scher [1,01 ■ More recently, Shlesinger wrote a review that 
contributed to further popularize the CTRW @ ; theoret- 
ical, numerical, and empirical studies on the CTRW have 
been discussed by Weiss , Metzler and Klafter [1, Q , 
and some authors of the present paper [l^, . 

In a CTRW, if X{t) denotes the position of a diffusing 
particle at time t, = X{ti)^ X{ti^i) denotes a random 



jump occurring at a random time ti, and = ti — ti^i 
is the waiting or sojourn or interarrival or duration time 
between two consecutive jumps, one has 



N{t) 



N(t) 
def \ ^ 

i=i 



(1) 



where to = 0, -'^(O) = and N{t) is a counting random 
process that gives the number of jumps up to time t. 
Throughout this paper, we assume that 

- the jumps ^i, i = 1,2,... are independent and 
identically distributed (iid) random vectors in R'^, 



1,2,... [Ill 



the waiting times r^, i 
variables in R_|_; 



the families (^i, i 
are independent. 



1,2,. 



1,2,... are iid random 
.) and (Ti, i = 1, 2, . . .) 
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The third assumption means that we consider a so-called 
uncoupled CTRW. The first two assumptions entail that 
the joint distribution of any pair {^i,Ti) does not depend 
on i. If, in the uncoupled case, the law of (fi, Ti) is given 
by a density function ip{^, t), the independence of and 
Ti means that it can be factorized in terms of the marginal 
probability densities for jumps A(^) and waiting times 
V'(r): v'(e,r) =A(C)V(t). 

Eq. Il]) means that a CTRW is a random sum of in- 
dependent random variables. The process of the jump 



2 



times 



Thus, if B = (0,t]: 



tn = ^ ■''t: to ~ Oi 



(2) 



1=1 



is a renewal point process. Therefore, a CTRW can 
be seen as a compound renewal process [13, [13, 
The existence of an uncoupled CTRW can be proved, 
based on the corresponding theorems of existence for 
renewal processes and discrete-time random walks [l5| . 
Cadlag (right-continuous with left limit) realizations of 
a CTRW can be easily and exactly generated by Monte 
Carlo simulation and plotted [l^l- This is illustrated in 
Fig. [TJ An uncoupled CTRW is Markovian if and only if 
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FIG. 1: Realization of a CTRW with exponentially dis- 
tributed waiting times {'yt — 1) and standard normally dis- 
tributed jumps = and a = 1). 

the waiting time distribution is exponential, i.e. iP{t) = 
exp(-T/7t)/7t [ll,[i3- An uncoupled CTRW belongs to 
the class of semi-Markov processes [13, [H, [H, [131, i-e. 

for any A C M'* and t > Q we have 

P{Sn e A,Tn < t \ So, . . . , Sn-l,Ti, . . . , r„_i) 

^P{SneA,Tn<t\Sn-l) (3) 

and, if we fix the position Sn-i = y of the diffusing 
particle at time tn-i, the probability on the right will 
be independent of n. In the generic coupled case, if the 
law of {£,mTn) is given by a density function (^^(Ct), we 
can use Sn ~ Sn—i + and rewrite this as 



ip{x — S'n-i, t) drdx. 



P{Sne A,T„<t\Sn-l) 

J A Jn 

(4) 

This can be shown as follows. Let Ia{x) denote the in- 
dicator function that yields 1 ii x G A and otherwise. 
Probabilities can be replaced with expectations writing 
P{x G A) = E[Ia{x)]. Moreover, one has IaIb = Iahb- 



P(5„eA,T„eB|5„_i) 

= S[/a(5„)/s(t„)|5„_i] 

= E[lAiS„-l+UlBiTr,)\Sn,l] 



f r lA{Sn-i+OlB{rM^,r)dTd^ 
Jw Jo 

[ [ lA{Sn-i+Ov{^,r)dTd^ (5) 

jRd J B 

t 

Ia{x)ip{x — Sn-i, t) drdx 
ip{x — Sn~i, t) drdx. 



A Jo 

MontroU and Weiss wrote Eq. (jlj as an integral equa- 
tion for the probability density px{x,t) of finding the 
particle in position x at time t in terms of the joint prob- 
ability density ip{£,, r) of the jumps ^ and waiting times 
r: 



Px{x,t)=5{x)-^{t)+ / ^{i,T)px{x~tt-T)dTdt 

JW Jo 

(6) 

where ^'(t) = 1 — Jg iP{t) dr is the complementary cu- 
mulative distribution function for the waiting times, also 
called survival function. This can be shown observing 
that 



P{X{t) e dx I X{Q) = 0) = px{x, t) dx 



(7) 



and 



P{X{t) e dx\X{t') 

= P{X{t - t') e dx I X{Q) = x ) 

= P{X{t-t')-x' C,dx\X{Q)^Q) (8) 

= Px{x — x', t — t') dx 

because the increments in time and space are iid and 
hence homogeneous. Moreover, from Eq. (|4]), 

P{Si e dx, Ti e dt\So = 0) = ip{x, t) dxdt. (9) 

The probability in Eq. ([7]) can be decomposed depending 
on the duration of the first jump ri with respect to t: 



P{X{t) e dx\X{Q) = 0) 

= P{X{t) (^dx,Ti>t\X{Q)^Q) 
+ P{X{t) C,dx,Ti<t\ X{Q) = 0). 

The part without a jump before t is given by 



(10) 



P{X{t) G dx, ri >< I X{Q) = 0) = P(ti > t)6{x) dx 

^ 5{x)'^{t)dx. (11) 
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The other part is given by 

Pix{t) e dx,Ti <t\x{o) = 0) 

= f f P{X(t) edx\X(t') 
Jr'' Jo 



Eq. ([gj 

Eq. lHJ 



X P{Si e dx',Ti G dt'lSo = 0) 

/ / P{X{t) e dx\X{t') ^ x')(p{x\t')dt'dx' 
Jm'' Jo 



Px{x — x\t — t')dx (p{x' , t') dt'dx' 



WL'^ JO 



JO 



t) pxix-tt- t) drdi 



dx. 



(12) 



Combining Eqs. ((IT|) and yields Eq. dH). Notice that 
the latter just gives a one-point probability density, which 
is not enough to characterize a stochastic process without 
further assumptions. 

Eq. ([6]) can be solved in the Fourier-Laplace domain, 



p{k,s) 



1 



l-V(s) 



(13) 



1 — (p{k, s) ^ 
where the Fourier and Laplace transforms are defined as 



+ 00 



/(x)e''^'^ dec, fceR, (14) 
f{t)e-'*dt, see. (15) 



The inverse transforms to the space-time domain are pos- 
sible in the uncoupled case, i.e. when (p(^,r) = X{£,)iP{t); 
this leads to a series expression written in terms of the 
probability P{N{t) = n) = pM{n,t) of the counting 
process N{t), and the 7i-fold convolution A*"(a;) of the 
marginal probability density of jumps A(^): 



px(a;,i) = ^Piv(r^,t)A*"(a;). 



(16) 



n=0 



The method using integral transforms is described in sev- 
eral papers, including the original one by MontroU and 
Weiss. However, Eq. ([T5|) can also be derived directly by 
probabilistic considerations. Indeed, Eq. ([T]) is a random 
sum of iid random variables. This means that any po- 
sition X can be reached at time i by a finite number n 
of jumps. The probability of reaching position x at time 
t in exactly n jumps is pAr(n, i)A*"(.T). Eq. pB]) follows 
given that these events are mutually exclusive. Note that 
Pjv(0, t)\*^{x) coincides with the singular term 5{x)'^{t), 
meaning that the distribution function for x has a jump 
at position x = of height 5'(t). 

A CTRW with exponential waiting times is called a 
compound Poisson process (CPP), as in this case 



pjv(n,i;7t) = exp(-t/7t) 



(17) 



A CPP is not only a Markov, but also a Levy pro- 
cess. This means that it has independent and time- 
homogeneous (stationary) increments. In the Levy case 
px(x,t), even px{x,l), fully characterizes the stochas- 
tic process defined by Eq. ^ [l^, [2l|, this is due 
to the infinite divisibility and the fact that the incre- 
ments are stationary and independent. For a normal 
CPP, i.e. a CPP with normally distributed jumps, the 
ri-fold convolution A*"(a;) of N{ii,a'^) can be evaluated 
as N{Tiii,n(j'^), leading to 



px{x,t;fi,a,jt) =exp(-t/7t) 



E 

Tl = 



exp 



{x — n^y 



(18) 



B. The CTRW in physics, insurance, finance, and 
economics 

Since the seminal paper by Montroll and Weiss [l|], 
there has been much scientific activity on the application 
of the CTRW to important physical problems. One line 
of research investigated anomalous relaxation related to 
power-law tails of the waiting time distribution as well as 
the as ymp totic behaviour of the CTRW for large times 
[1, [11,111 m, SHI- As mentioned above, Klafter and 
Metzler have extensively reviewed these and subsequent 
studies [1, Q. Furthermore, in their book, ben-Avraham 
and Havlin have discussed the applications to physical 
chemistry [28| . Here, it is worth mentioning the recent 
work on the relation between the CTRW and fractional 
diffusion that can be traced to papers by Balakrishnan 
and Hilfer [2 9|. [S Ol and has been thoroughly discussed in 
Refs. [l0l . l3ll . l32| . Some specific applications include, e.g., 
plasmas [33| and biopolymers [3J, |35| . 

The CTRW has been applied also in insurance, finance, 
and economics. Even if well-known in the field of econo- 
physics [s^ . [36| , these applications deserve a short sum- 
mary. 

In ruin theory for insurance companies, the jumps 
are interpreted as claims and they are positive random 
variables; ti is the instant at which the i-th claim is paid 

na. 

In mathematical finance, if PA{i) is the price of an 
asset at time t and -Pa(O) is the price of the same as- 
set at a previous reference time Iq = 0, then X[t) = 
\og{PA{t) / Pa{^)) represents the log-return (or log-price) 
at time t. In regulated markets using a continuous 
double- auction trading mechanism, such as stock mar- 
kets, prices vary at random times ti, when a trade takes 
place, and = X[U) - X{U_i) = \og{PA{U) / Pa{U-i)) 
is the tick-by-tick log-return, whereas = ti — i^-i is the 
intertrade duration; for more details, see [3^. [36l. [sst and 
references contained therein. 

In the theory of economic growth, represents a 
growth shock, which can actually be both positive and 
negative, X{t) is the logarithm of a firm's size or of an 



4 



individual's wealth, and is the time interval between 
two consecutive growth shocks; see [s^] and references 
therein. 



C. Motivation for the study of stochastic integrals 
driven by a CTRW and link with fractional calculus 

Given the wide range of applications of the CTRW 
ovcrviewcd in the previous subsection, it is relevant to 
study diffusive stochastic differential equations whose 
driving noise is defined in terms of a CTRW: 



dZ = a{Z, t)dt + b{Z, t)dX. 



(19) 



Here Z{X, t) is the unknown random function, a(Z, t) 
and t) are known functions of Z and time t, and dX 
represents the CTRW 'measure' with respect to which 
stochastic integrals are defined. In order to give a rig- 
orous meaning to such an expression, some constraints 
on the properties of the CTRW are necessary. In a re- 
cent paper, the theory has been discussed for stochastic 
integration on a time-homogeneous (stationary) CTRW 
— i.e., the already mentioned CPPs [sl]. Although the 
theory reported there was already well known by math- 
ematicians and has been used in finance for option pric- 
ing since 1976 [40| . that paper contains useful material 
and is written in a way that is clear and appealing for 
physicists. Here, inspired by Ref . [s^ , the theory will be 
further discussed and developed. 

Consider a CTRW X{t) whose jumps in space are 
distributed according to the symmetric Levy a-stable 
law, a e (0, 2], whose density can be expressed as a series 
or, more conveniently, as the inverse Fourier transform of 
its characteristic function: 



ia(C;7x)=^,7Mcxp(-|7.fcr)] (0- 



(20) 



For a = 2 this corresponds to a Gaussian with standard 
deviation a = v^Tx- Let the waiting times of the 
CTRW have the probability density 



^l3{T\lt) 



dT 



-Er, 



(21) 



where Ep^z), (i G (0, 1], is the one-parameter Mittag- 
LefSer function [ill. 142. [43| : 



z e 



(22) 



For a real argument z = t e R and /3 = 1 this corresponds 
to an exponential function. When < 1, Epl—t^) is ap- 
proximated for small values of i by a stretched exponen- 
tial decay (WeibuU function), exp {-t'^/r{l + (3)), and 
for large values of i by a power law, t~'^/r(l — /?). 

In the diffusive limit for X{t), when the scale pa- 
rameters "fx of the jumps and 74 of the waiting times 
vanish satisfying the scaling relation j"/jt = D, if 



in Eq. ([19)) a = and 6=1 the probability density 
Pz{z, t) = px{x, t; 7a;, 7t) converges to the solution of the 
space-time fractional diffusion equation (FDE) [3, 



gl3 ga 

—ux{x,t-D)=D-^ux{x,t-D) 



(23) 



ux{x,Q'^-,D) = 5{x), x£ 



t g 



The space-fractional derivative of order a E (0, 2] is de- 
fined according to Ricsz: 



/(^) = 



(24) 



The time-fractional derivative of order f3 <E (0, 1] is de- 
fined in the sense of Caputo 



dP_ 

dtl^- 



(25) 



The FDE is a generalization of the standard diffusion 
equation, that results for a = 2 and /3 = 1; in this case 
the solution ux{x,t] D) of the Cauchy problem given 
by Eq. (|^ is the one-point probability density of the 
Bachelicr- Wiener process or Brownian motion i?(t). 



ux{x,t; D) 



1 



exp 



x^ \ 
4Dt) ' 



(26) 



and X{t) is the NCPP introduced at the end of SecIS 
The general solution of the FDE was worked out in the 
Fourier-Laplace domain: 



uxik,s) 



D\k\' 



(27) 



Because 



c: 



^0-1 



Dlkl" 



it) = Ep{-D\krt^) (28) 



defining k = kt^^" and the time-independent Green func- 
tion 

Go^A^; D) = [Ep{-D\Kn] (0, (29) 

can be expressed in 



(30) 



the solution of the FDE, Eq 
the space-time domain as 

ux{x, t; D) = t-P/^ G^Axt-P/^- D). 



These results arc a consequence of a generalized central 
limit theorem for sequences of random variables [sij . A 
simpler derivation can be found in Ref. 3^. For com- 
putational details see Sec. IIIII and Ref. 10]. If a{x,t) 
and b{x,t) are not constant, a fractional Fokker-Planck 
equation ior ux (x,t; D ) has been proposed in the diffu- 
sive limit [sl. [461. [47[. [Isl. l49l. Isot starting from a generalized 
master equation (^jor a CTRW Hi]. For the NCPP this 
reduces to the standard Fokker-Planck equation [si], H^] ■ 
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Without taking the diffusive hmit, and if a = and 
b ~ 1, the time evolution of the probability density 
Px{x, t) is given by the Montroll- Weiss integral equation 
([5]). The uncoupled case of the latter can be presented 
alternatively in an integro-differential form (53j , 

/■* d 

J ^{t - T)—px{x,T)dT 

/ + 00 
Xix~OPxiC,t)d(, (31) 
-oo 

that can be interpreted as a time evolution equation of 
Fokker-Planck type. It involves the time derivative of 
px{x,t) and an auxiliary function <i>(t) defined through 
its Laplace transform as $(s) ~ \D'(s)/V'(s), so that 
'^[t) = Jp $(<: — t)iP{t) dr. This approach has been gen- 
eralized studying scores of possible kinetic equations for 
non-Markovian processes [5J| . What follows in the next 
sections is valid without necessarily taking the diffusive 
limit. Nevertheless, the latter is important because it 
motivates our particular choice for the marginal distribu- 
tions of jumps and waiting times, and because it provides 
analytic expressions that can be compared to our Monte 
Carlo results as shown in Sec. IIIIl 

II. STOCHASTIC INTEGRALS 

In Ref. [3§], the stochastic integral is never explicitly 
defined. However, starting from the fact that sample 
paths of a CTRW can be represented by step functions, 
it is possible to give an explicit formula. 

A. Definitions 

Some heuristic manipulations are useful for the defini- 
tion of the stochastic integral 

J{t)= f Y{s)dX{s), (32) 
Jo 

where X{t) and Y{t) arc synchronous CTRWs, i.e. their 
jumps happen at the same times U, i = l,...,N{t). 
Though an interesting case is often Y{t) = G{X{t)) with 
a suitable function G{X), the jumps of Y{t) and X{t) at 
t = ti may be independent as well. Eq. ([T]) defining X{t) 
can be written in terms of the right-continuous variant 
of Heaviside's step function 9{t), which is for t < and 
1 for t > 0: 

Nit) 

X{t) = mt - U). (33) 

1=1 

Using the fact that the 'derivative' of Heaviside's 9 func- 
tion Oit — ti) is Dirac's 5 function 5{t — ti), one can write 

Nit) 

dX{t)^Y,i,5{t-t,)dt, (34) 

i=l 



which means that AX(ti) =^ X{ti) — X{t^ ) = with 
X{t^) = lim^^j- X{s) lim,^t,,^<t, X{s). Note that 
5{t) is not a proper function, but rather a distribution in 
the sense of Sobolev and Schwartz [sll . Writing Eq. ([M]) 
with t~ in place of ti, inserting it into Eq. (j32p . and using 
the properties of Dirac's 5 function, we get the exact 
expression (no limit needed: recall that the number of 
jumps N{t) between and i is a random finite integer) 

l-t Nit) 

m / y(s-)dx(s) = ^r(t-)6 

Nit) 

= Y.Yitr){X{t,)-X{tr)). (35) 

i=l 

The choice Y{s^) for the integrand makes I{t) a mar- 
tingale if X{t) is a martingale, as will be explained be- 
low. This naive definition works nicely if the driving 
noise is a step function with jump times ti and jumps 
£_i = X{ti) — X{t~); if Y{t) and X{t) jump at the same 
time we even have Y{t~) = Y{ti^i). As soon as one 
wants to go beyond this situation, measurability and con- 
vergence become an issue. This observation prompted 
K. Ito to use martingale convergence theorems to tackle 
the convergence for a large class of integrators [s^ . To do 
so wc must make sure that I{t) is a martingale whenever 
X{t) is. For this we assume that Y{t) is adapted i.e. mea- 
surable with respect to the natural filtration generated by 
the driving noise: Tt = cr{X{s) : s < t). Therefore the 
integrand Y{t^) in Eq. ([55]) becomes statistically inde- 
pendent of the increment ^i = X(ti) — X[t~) and we end 
up with a stochastic integral I[t) that is a martingale; 
see the next section for details. The fact that we evalu- 
ate Y{t) at the left end-point t^ of the 'infinitesimal in- 
terval' [t^,ti] makes the integrand non-anticipating and 
adapted, i.e. independent of the increment. This can be 
seen as a causality requirement: one does not want Y{t) 
to anticipate the future behavior of ^(t) [S^l- An elemen- 
tary introduction to the concept of a non-anticipating 
function can be found in Ref. (ssj . Any adapted pro- 
cess with right-continuous (or left-continuous) paths is 
progressively measurable. 

In Eq. (|35p wc might equally well choose to evaluate 
Y{t) in the right end-point ti of the infinitesimal interval 
[t^ , ti] , corresponding to the right-continuous variant of 
Heaviside's 9 function in Eq. ([221), or in any intermedi- 
ate point if. This means, however, loosing the martin- 
gale property of the stochastic integral. The effects on 
the formulae for such a choice can be nicely described 
for random step functions Y{t) and X{t) jumping at the 
same times ti,t2, ■ ■ ■ ,tn- Write 

l-t Nit) 

Mt) - / r(s,)dx(s) = ^y(if)^, (36) 

•^0 ,= 1 

N(t) 

= 5] [(1 - d)Yit-) + ^Yit,)] [Xit,) - Xit-)] 
1=1 
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for a parameter G [0, 1] that interpolates linearly be- 
tween Y{t~) = Y{ti_i) and Y{ti), resulting in a con- 
tinuous class of stochastic integrals. The choice -d = 
gives the Ito integral Jo(i) = I{t). For any value of 
■& the integral is a right-continuous function with jump 
AJ^{U) = MU) - Mt-) = Y{tf )AXiU). 
Eq. ([55)) can be rearranged to 

(t) = J,/2 (f) + (^^ - [X, Y] it) , (37) 

where 

Nit) 

[X,Ym = Y.IX{U) - XitiWiU) - Yit-)] (38) 
1=1 

is the covariation or cross variation of X{s) and Y{s) 
for s e [0, i]. When Y{s) = X{s), the quadratic 
variation [X, X] (t) is denoted simply as [X] (t) . Thus 
each member of the family of stochastic integrals with 
1? £ [0, 1] can be obtained adding a compensator to the 
Stratonovich integral Ji/2{t) = S{t). The latter corre- 
sponds to the symmetric variant of Heavisidc's step func- 
tion, 9{t) = (sgnt -I- l)/2, and is particularly appcahng 
because it can be computed according to the usual rules 
of calculus. However, the Ito integral has the advantage 
of being a martingale, as proved in the next subsection. 
The distinction between integrals with different values of 
•& disappears in the continuous limit for processes with 
finite variation, e.g. continuously differentiable functions, 
because this implies that their quadratic variation is zero 
[5^. Unless stated otherwise / Y{s)dX{s) denotes the 
Ito integral, while the Stratonovich integral is often indi- 
cated as / Y{s) o dX{s). 

B. Martingale property of the Ito integral 

Although it is easy to simulate directly the stochas- 
tic process defined in Eq. (P5|) — see the next section 
for numerical examples — it is not so easy to derive its 
properties. Each term in the sum depends on the previ- 
ous ones and the nice properties of convolutions are not 
helpful here. However, using the martingale transform 
theorem, it is possible to obtain conditions under which 
I(t) is a martingale. 

In order to define martingales, we need a filtered prob- 
ability space {n,J^,{Tt)t>o,P), where {Tt)t>o is a fil- 
tration — i.e., an increasing family of sub tr-algebras — 
representing the information available up to time t. A 
martingale is a stochastic process X{t) for which the ex- 
pected value E[|X(t)|] exists for t > and the conditional 
expectation E[X{t) \ J^s] is X{s) for alH > s [H, Hi, HO] ■ 

Let us consider the natural filtration, that is the a- 
algebra generated by the CTRW itself: Tt = (j{x{s) : 

s < t) ^ cr(Ci, . . . ,^fc;Ti, . . . ,rfc : k < N{t)) =^ GN{t}- 
Then X(t) is a martingale with respect to J-t if and only 
if the mean of the jumps is zero. Denote by {ti,S,i) 



the time and height of the finitely many jumps i = N(s) + 
1, . . . , N{t) occurring between s and t > s. Then 

N{t) 

E[Xit)\T,]^X{s)+ lEfel.f's]. (39) 

i=JV(s) + l 

Using the semi-Markov property, Eq. ([3]), we get for i > 

Nis) 

E[e. I = E[e. I GNis)] - m I ^Nis)] = m] - o, (4o) 

thanks to the independence of and Ci, . . ■ ,^jv(s)- 
Eq. ((Ml) becomes 

E[X{t)\J^s]^X{s), (41) 

which shows that {X{t))t>o is indeed a martingale with 
respect to its natural filtration. 

Note that our argument is valid for a general uncou- 
pled CTRW. We do not need the independence of the 
increments X{t + At) — X{t) of the process X{t) for non- 
overlapping intervals. Of course, if we have independent 
increments, i.e. a compound Poisson process X{t), the 
proof becomes easier. 

Let us now investigate the integral defined in Eq. ([55]) 
for a martingale CTRW X{t). If there is an arbitrary 
but finite number of jumps between s and t > s, one has 

N{t) 

E[Iit)\T,]^I{s)+ nY{t~)^^\gNis)]; (42) 

i=N{s) + l 

now, one observes that = X{ti) — X{ti-i) and that the 
random sum in Eq. (|42p becomes 

N{t) 

J2 E[r(tr)e,|5^(,)] 

i=N{s) + l 

N{t) 

= Y E[F(t-)(X(tO-X(t,_i))|5^.(,)]. (43) 

If Y{t) is measurable with respect to !Ft = GN(t), then 
Y{t~) is -(-measurable. Since N{t~) = N{ti-i), 

this means that Y(t^) is fJjY(ti_i) = ^i-i-measurable; 
this is to say that Y{t~) is predictable for the filtration 
Qi, i.e. the value of Y{t~) is known at time t^-i. When- 
ever for each i the expression Y{t~){X{ti) — X{ti^i)) 
has a finite absolute mean — e.g., if the process Y{t^) is 
bounded — we have 

E[y(^-)(x(t,)-x(^,_l))|e^^(,)] 
= E[E[Y{t;){x{t,) - x(t,_i)) I I g^^,)] 

= E[y(t-) E[(X(t,) - X(t,_i)) I I 

In the above calculation we have used the fact that Gn^s) 
is contained in Qi-i as (z — 1) > N{s), along with the 
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tower property and the fact that we can pul l out what is 
known from the conditional expectation [60|. Since X{t) 
is a martingale, we have M[X{ti) \ ^ti-i] — X(ti-.i) which 
means that 

E[Y{t-){Xit,) - X(t,_i)) I Gms)] = 0. (45) 

Consequently, each term in the random sum vanishes and 
E[/(t) I !Fs] = -f(s)- Summing up, if X{t) is a martingale 
with respect to J-'t and if the integrand is bounded and 
predictable, one has that I{t) is also a martingale with 
respect to Tf 



III. SIMULATION 

In the previous section we have explicitly defined and 
rigorously characterized a martingale stochastic integral 
driven by an uncoupled CTRW and given in Eq. ([55]) . as 
well as a more general class of stochastic integrals given 
by Eq. ((5^ . A useful property of these equations is that 
they can be easily implemented by means of Monte Carlo 
simulation, as will be shown here for the case Y{t) = 
X{t). The theory of Sec. |TT] is the basis for the Monte 
Carlo solution of stochastic differential equations driven 
by CTRWs and discussed above in Sec. II CI 

The marginal distributions of jumps and waiting times 
presented in Sec. II Cl are apparently demanding, but they 
can be sampled easily using one-line transformation for- 
mulas [lol . [6ll . [6^ . A random number ^ drawn from the 
symmetric Levy a-stable probability density, Eq. (PO)) . 
can be obtained from two independent uniform random 
numbers U,V ^ (0, 1) through a transformation due to 
Chambers, Mallows and Stuck [il,!!!, 



— log U cos $ 

C0S((1 - Q!)$) 



sin(Q!$) 
cos $ 



(46) 



where $ = irlV - 1/2). For a = 2 Eq. (gHl) reduces to 
^ = 2jxV~ logf/ sin<I>, i.e. the Box-MuUer method for 
Gaussian deviates with standard deviation a = \/2jx ■ A 
random number t drawn from the one-parameter Mittag- 
Leffler probability density, Eq. (|2ip . can similarly be 
obtained from two independent uniform random num- 
bers U,V € (0,1) through a transformation proposed by 
Kozubowski and Rachcv [65,, £61] : 

For (3—1 Eq. (|47p reduces to the transformation formula 
for the exponential distribution, r = —74 log U . 

Now, as outlined above, the Monte Carlo simulation of 
an uncoupled CTRW is straightforward. To compute the 
value X{t), generate a sequence of N{t) + 1 iid waiting 
times Ti until their sum is greater than t. Discard the 
last waiting time and generate N{t) iid jumps ^i. Their 
sum is the desired value of X{t). Based on Eqs. ^ and 
(12); this algorithm was used to generate Fig. [TJ This 



procedure is also the basis to compute I{t) according to 
Eq. ([55)) . or more in general J^{t) according to Eq. 
and the covariation [X, Y\{t) according to Eq. ([38|) . Each 
jump ^, is multiplied by Y{t-), {l-'d)Y{t-) + W{ti), or 
Y{ti)~Y{t~), and the results of these multiplications are 
summed to obtain respectively I{t), J§{t) and [X, y](t). 
C-f-f- code for the case Y{t) = X{t) can be found in the 
appendix. 

Figs. [2] and [3] show histograms from 1 million Monte 
Carlo realizations ofX{t), I{t) = J*X{s)dX{s), S{t) = 

JgX{s) o dX{s) and [X]{t), where i = 1 and X{t) is a 
symmetric CTRW with jump and time scale parameters 
linked by the relation j^/j^=D = l. Thus the integrals 
in Figs. [2] and [3] give the Monte Carlo solution for t — 1 
of the stochastic differential equation dZ = XdX with 
initial condition Z(0) = 0. Since the Ito integral is a mar- 
tingale starting at zero, its mean is zero. This is not true 
for the Stratonovich integral. The probability density of 
the Stratonovich integral S(t) = X'^{t)/2 can be worked 
out from the density of the stochastic process X{t) by 
the transformation p5(s, t) = J2iPxixi{s),t)\dxi{s)/ds\, 
where the sum is over all Xi that yield the same s. For 
s ~ this is X1.2 — ±-\/2s and thus 



psis,t) = 2pxiV2^,t)/V2^, s>0. 



(48) 



In the diffusive limit the NCPP X{t) approximates the 
Bachelicr- Wiener process B{t) [s^, and thus the proba- 
bility density of the process X{t) approximates the den- 
sity of B{t), Eq. ((26)) . The analytic probability den- 
sity for the Stratonovich integral S in the diffusive limit 
can be obtained inserting the probability density of the 
Bachelier- Wiener process into the transformation for- 
mula given by Eq. (|48)) . yielding 



Ps{s,t;D) 



, cxp f--^-) , s>Q. (49) 

V2^^ V 2DtJ ^ ' 

According to Eq. §7} here /(t) = S{t)~[X]{t)/2; if the 
dependence of S and [X] is small, the probability density 
of the Ito integral is approximated by the convolution of 
the probability density of the Stratonovich integral with 
that of the quadratic variation mirrored around zero and 
scaled to half its width: 



Pi{x,t) 



+ 00 



ps{x + 2x',t)p^x]{-'2x',t)dx'. (50) 



For all choices of a and j3 the agreement between the 
analytic expressions for X{t) and S{t) in the diffusive 
limit and the empirical results from Monte Carlo simu- 
lation of the CTRWs is fair already for the largest value 
74 = 0.1: the curves cannot be distinguished by eye 
at the scale of our plots. Therefore we did not evalu- 
ate the analytic probability density for X{t), Eq. p^ . 
available for the particular case of a NCPP only, i.e. 
the left column of Fig. [21 Instead the quadratic varia- 
tion [X] (t) and consequently the Ito integral tend visibly 
more slowly to their diffusive limits. For a NCPP the 
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FIG. 2: (Color online) Convergence of the empirical probability densities p from 1 million Monte Carlo runs (points) to the 
analytic probability densities u (lines) in the diffusive limit for a CTRW X{t), its Stratonovich integral S{t), its Ito integral 
I{t), and its quadratic variation [X](t), with t — 1 and different choices of the index parameters a, f3 and of the scale parameters 
7a;, 7t, where "f^ /'Jt = D = 1. E[Af(f)] is the average number of jumps per run. 
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diffusive limit of [X]{t) is [B]{t) = 2Dt. In this limit 
I[t) = S{t) - Dt = B^{t)/2 - Dt, corresponding to the 
well-known result that the probability density of the Ito 
integral is equal to the density of the Stratonovich inte- 
gral shifted by —Dt, i.e. pi{x, t) = ps{x + Dt, t). Though 
the quadratic variation of the NCPP is appreciably dif- 
ferent from its limit 5{x ~ 2Dt), where Dt = 1, for any 
non infinitesimal value of 74 as shown in the left column 
of Fig. [21 for 7t = 0.01 there is a good agreement between 
the Ito integrals from Monte Carlo and from Eq. ([50]) . 

The density of the quadratic variation for a CTRW 
can be obtained from the density of squared jumps, 
\(^i{x), that results from a transformation of the den- 
sity of jumps, \^{y/x), similar to the one that leads from 
px{x,t) to ps{x,t), Eq. (|48p . except for a factor 2: 



(51) 



Inserting this equation into the solution of the MontroU- 
Weiss equation in the space-time domain, Eq. (|16p , gives 



Pix]{x,t;-f^,jt) = ^PAr(n,t;7t)A^?(x;7^), (52) 



n=0 



where a; > 0. Unfortunately even for an NCPP the n-fold 
convolution cannot be computed as easily as for px{x, t) 
in Eq. p8| . However, the characteristic function of the 
quadratic variation can be written as 

00 

P[x]{k,t]-f.^,-1t) = ^PAr(n,i;7t)A^2(A:;7a:). (53) 

n=0 

In order to consider non-exponential waiting times with 
power-law tails and infinite first moment, for the sake of 
simplicity let us assume that piq (n, t) is the distribution 
of the Mittag-Leffler counting process [3l| , 



Pn 



where 



E 



(») 

/3 



(54) 



(55) 



This choice is more general than it seems, as the Mittag- 
LefSer distribution for waiting times is an attractor for 
the thinning procedure used to obtain the diffusive limit 
(67| . Using the Mittag-Leffler distribution from the be- 
ginning simplifies the derivation of this limit. Then 
Eq. (IS2D becomes [32| 



A^2(fc;7,))). (56) 



As the jumps ^ follow a Levy a-stable distribution, for 
X ^ 00, X^2{x;^x) ~ ix/"fx)~"^^~^ , and the sum of 
converges to the positive stable distribution with index 
a/2, whose characteristic function is 

A42(fc;7,) = i;+/2(fc;7.) ^ exp {i-^^lkr/^). (57) 



The scale parameter is the same as in the Levy sta- 
ble distribution, Eq. (PP]) . Inserting this distribution in 
Eq. (|56p . the diffusive limit yields the following charac- 
teristic function for the quadratic variation: 



{k,t]D) ^ Ep{- D{-ik)"/H^) . 



(58) 



Now we can proceed in a similar fashion as for the so- 
lution of the FDE, Eqs. (pSHgO)) . Defining k = kt'^f^/" 
and 

M^A^;D) = T-'[Ep{ - Di-zKr/')]{0, (59) 

where ^ > 0, we obtain the quadratic variation for the 
diffusive limit in the space-time domain. 



U[x]{x, t- D) = t-2/3/« M„,^(xi-2/3/"; i^). 



(60) 



When a — 2, M2./3(^) coincides with the right half of 
the Mainardi- Wright function [gl], which is also called 
M-function of Wright type because its shape recalls a 
capital M centered in the origin. When a = 2 and (3=1 
(standard diffusion case), a delta function u\^x] {^i D) = 
6{x — 2Dt) is recovered, corresponding to the quadratic 
variation of the Bachelier- Wiener process, [X] {t) = 2Dt. 
The plots in Figs. [2] and [3] display quadratic variations 
both from Monte Carlo and from Eq. (|SD1) . The conver- 
gence of the quadratic variation in the diffusive limit can 
be used to prove that the integrals of X(t) as defined in 
Sec. lU converge. 



IV. CONCLUSIONS AND OUTLOOK 

This paper is based on the definition, given in Eq. (|36p . 
of a class of stochastic integrals J^{t) driven by a CTRW 
X{t). For = this results in the Ito integral I{t), 
Eq. ((35)l . for d = 1/2 in the Stratonovich integral. If the 
process X{t) that defines the measure used in Eq. ([551) 
a martingale with respect to its natural filtration, then 
I{t) is a martingale too; this is a consequence of the mar- 
tingale transform theorem. It turns out that an uncou- 
pled CTRW with zero-mean jumps is a martingale. The 
stochastic integration theory developed here is more gen- 
eral than the one sketched in Ref. [39[ , as it can be applied 
also to a CTRW that is neither Markovian nor Levy. In 
fact, exponential waiting times are not needed to prove 
that I{t) is a martingale if X{t) is a martingale. 

The theory presented in Sec. [Ill lies at the foundation 
of the Monte Carlo method for integrating stochastic 
differential equations driven by CTRWs. As explained 
in Sec. [H these results are relevant for applications in 
physics and economics as well as in all those fields like in- 
surance and finance where martingale methods can help 
in the quantitative evaluation of risk. Eq. p6p is a conve- 
nient basis for the Monte Carlo calculation of stochastic 
integrals. This is shown in Sec. IIIIl where Monte Carlo 
realizations of CTRWs are used to effectively approxi- 
mate the Ito and Stratonovich integrals driven by the 
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FIG. 3: (Color online) Convergence of the empirical probability densities p from 1 million Monte Carlo runs (points) to the 
analytic probability densities u (lines) in the diffusive limit for a CTRW X{t), its Stratonovich integral S{t), its Ito integral 
I{t), and its quadratic variation [X](t), with t — 1 and different choices of the index parameters a, f3 and of the scale parameters 
7a;, 7t, where "f^ /'Jt = D = 1. ¥,[N{t)] is the average number of jumps per run. 



11 



Bachclier- Wiener process and, more generally, by the so- 
lution of the space-time fractional diffusion equation. 

We believe that up-to-date mathematical methods 
from probability theory and stochastic calculus are bene- 
ficial to the study of the CTRW and of other random pro- 
cesses useful in statistical physics. We fear that progress 
will be slower or impossible if these methods are ignored 
by physicists. 

Future work will deal with Monte Carlo simulations 
for coupled CTRWs where jumps and waiting times obey 
fat-tailed distributions [gl, [T^I • There will also be a dis- 
cussion of convergence based on the results collected in 
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jumps = 

// Loop over runs 

for (run = 1; run <= runs; run++) { 

// Initialize and increment t, x, etc. 

t = 0, X = 0, qvar = 0, ito = 0, str = 0, 

tau = random. t(); // Eq. (47) 

while (t + tau < t_max) -[ 

t += tau; // time t 

xi = random. x(); // Eq. (46) 

qvar += xi*xi; // [X(t)] 

ito += x*xi; // I(t) 

str += (x+xi/2)*xi; // S(t) 

X += xi; // X(t) 

tau = random. t(); // Eq. (47) 

jumps++; // N(t) 

} 

// Update histograms at the end of each run 

hisx.add(x) ; // X(t) 

hisq. add (qvar ) ; // [X(t)] 

hisi.add(ito) ; // I(t) 

hiss.add(str) ; // S(t) 



Appendix 

Below are salient lines from the central loop of our 
C++ program for the Monte Carlo calculation of a 
CTRW X{t), its quadratic variation its Ito inte- 

gral I{t) and its Stratonovich integral S{t) as described 
in Sec. IIIII and shown in Figs. 



CPU times grow linearly with the number of jumps 
N{t) and take 1-3 fisec per jump depending on a and /3 
on a 2.2 GHz AMD Athlon 64 X2 "Toledo" Dual-Core 
processor with Fedora Core 7 Linux, using the Ran uni- 
form random number generator [t^ and the GNU C++ 
compiler (g++) version 4.1.2 with the -03 -static opti- 
mization options. 
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